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Using kinetic Monte Carlo simulations and a bond- 
counting ansatz, thermal stability and diffusion of an adatom 
island on a crystal surface are studied. At low temperatures, 
the diffusion constant D is found to decrease for a wide range 
of island sizes like D oc N~ a , where a is close to one, N be- 
ing the number of adatoms in the cluster. By heating up the 
surface, the system undergoes a phase transition above which 
the island disappears. Characteristics of that transition are 
discussed. 

PACS: 68.35.Fx, 82.20.Wt, 36.40.Sx 



I. INTRODUCTION 

Stability and dynamics of adatom as well as vacancy 
islands of monoatomic height on crystal surfaces have 
been studied intensively during the last years, both ex- 
perimentally and theoretically [1-12]. In particular, the 
equilibrium island size at low temperatures [10] as well 
as the decay of clusters in the presence of surface steps 
or larger islands ('Ostwald ripening') [9,11-12] have been 
analysed. Different microscopic mechanisms have been 
discussed to explain the size dependence of the diffusion 
constant D characterising the motion of the center of 
mass of an equilibrated island. Typically, D decreases 
with the number of adatoms N in the cluster in the form 
of a power-law, D oc N~ a , with a depending on the 
mechanism driving the island motion [1-8] . 

In this article, we consider adatom islands on a (100) 
surface of a simple cubic crystal with energy barriers for 
jumps of the adatoms to neighbouring surface sites given 
by isotropic nearest-neighbour bond energies. Equilib- 
rium and dynamic properties are computed by using ki- 
netic Monte Carlo techniques [13]. Specificly, for a single 
island in equilibrium with the 'gas' of adatoms on the 
surrounding terrace a phase transition is observed at a 
critical temperature T c , above which the island tends to 
disappear. This aspect seems to have been overlooked 
in previous investigations. A study on the equilibrium 
cluster size in a related Ising model had been performed 
already several years ago [14], but it had been limited 
to low temperatures. Quite recently, the thermal dis- 
integration of a cluster was discussed [15], however, for 
systems without conservation of the number of particles. 
In addition, we monitor the motion of the equilibrated 
adatom island well below T c , to estimate the value of a 
in the nearest neighbour isotropic bond-counting case. 



The article is organized accordingly. In the next Sec- 
tion, we introduce model and method. We then present 
results on the characteristics of the phase transition at 
which the island disappears, followed by a discussion of 
the diffusive motion of the cluster at low temperatures. 
We conclude with a short summary. 



II. MODEL AND METHOD 

Adatoms on a square lattice may be constrained to 
a single layer, with the occupation variable rii at surface 
site i being either 1 or 0. To jump to an empty neighbour- 
ing site, the adatom has to overcome an energy barrier 
E a . Using a bond-counting ansatz, that energy may be 
written in the form 



E a = E + nE b 



(1) 



where Eq is the activation energy for free diffusion of the 
adatom on a locally perfect surface, n is the number of 
occupied neighbouring sites, n— 0,1,2, or 3, and E b is the 
bond energy. Of course, the ansatz is not expected to give 
a realistic description of a specific material. However, it 
is useful in identifying generic properties of islands, as 
will be discussed below. 

The energy barrier (1) corresponds to the Hamiltonian 
(apart from a constant) 



n = s b 



TliUj 



(2) 



where the sum runs over neighbouring sites i and j. It 
is interesting to note that the Hamiltonian may be tran- 
scribed to a nearest-neighbour Ising model [11], 



Hx = -E b /A ^SiSj 



constant 



(3) 



with the Ising spin Si = ±1 being related to the occu- 
pation variable m by Si = 2ni — 1. The conservation of 
the number of adatoms on the surface corresponds to the 
conservation of the magnetization in the Ising model. 

We considered square surfaces with L x L sites and 
M x M adatoms, i.e. with a coverage 8 = M 2 /L 2 . Usu- 
ally, full periodic boundary conditions were used, to re- 
duce boundary effects. However, to compare with possi- 
ble experiments, we also applied free boundary conditions 
as encountered for a terrace of L 2 sites with large reflect- 
ing energy barriers at the descending steps bordering the 
terrace (large Schwoebel-Ehrlich barriers). 

In the simulations, L ranged from 50 to 1000, and M 
from 5 to 50, with the coverage 9 varying between 4 x 10~ 4 
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and 4 x 10~ 2 . Obviously, at those coverages the adatoms 
form a compact island at low temperatures due to the 
attractive interactions, see (1). A typical configuration 
is shown in Fig. 1. 

The kinetic Monte Carlo simulations were performed 
in the standard way [13,16], based on jump probabilities 
for the adatoms to neighbouring sites oc exp(— Ea/ksT), 
with ks being the Boltzmann constant and T the tem- 
perature. The time may be measured in units of trial 
jumps per adatom (MCA). Other commonly used time 
scales, invoking, e.g., a microscopic attempt frequency v, 
are linearly related to that unit. At low temperatures, 
the efficient algorithm of Bortz, Kalos, and Lebowitz [17], 
BKL, was implemented. 

To study stability and dynamics of the island, several 
quantities were computed. Among others, we recorded 
the distribution of clusters (as usual, a s-cluster is formed 
by s adatoms connected by occupied neighbouring sites), 
especially the fraction of adatoms in the largest cluster, 
i.e. the 'reduced island size', 

Umax — N max /M (4) 

where N max is the number of adatoms in the largest clus- 
ter. Furthermore, the fluctuations of that island size, the 
density of adatoms, the energy E, see (2), and the fluc- 
tuations of the energy were monitored. To analyse the 
motion of the island, we determined the time evolution 
of the position of its center of mass, see Section IV. 

III. THERMAL STABILITY 

In the ground state, T = 0, the M 2 adatoms form 
a compact square cluster, at sufficiently low coverage 9 
(otherwise, a stripe of adatoms will minimize the energy). 
At non-zero temperatures, adatoms may detach from the 
island, leading to a dynamic equilibrium of the rounded 
cluster and the 'gas' of adatoms on the terrace, see Fig. 
1. Of course, the size of the island is expected to shrink 
as the temperature increases, as demonstrated in Figs. 2 
and 3. 

For reasons of simplicity, we assume equal energy barri- 
ers in the bond-counting ansatz, (1), i.e. E = E b . Start- 
ing the simulations with a square cluster of M 2 adatoms, 
the time evolution of the island size, Tiraaxit) , at various 
temperatures ksT/Eb, is shown in Fig. 2. Obviously, 
equilibration may be rather slow, as observed before in 
the related Ising model [14] . Discarding the initial relax- 
ation, the thermal equilibrium value for the reduced is- 
land size, n m (T) = (n max ), is obtained by averaging over 
the subsequent, possibly strongly fluctuating, see Fig. 2, 
simulational data. The resulting temperature dependent 
island size is depicted in Fig. 3, at coverage 9 = 0.04 
with M= 25 and 50 (hence, L= 125 and 250). 

The drastic decrease of n m (T), both for periodic and 
free boundary conditions, in a narrow range of temper- 
atures suggests a phase transition at the critical tem- 
perature T c in the thermodynamic limit, L, M — ► oo, 



at constant coverage 9 = M 2 /L 2 . The reduced island 
size n m (T) may then be interpreted as the order pa- 
rameter, vanishing at T > T c . In the high-temperature 
phase, the island disappears in the gas of adatoms. In- 
deed, finite-size analyses of the Monte Carlo data, at 
fixed coverages, allow to locate the phase transition: For 
example, at 6 = 0.01, the turning point, T m (M), of 
the temperature dependent island size n m (T), see Fig. 
3, is found to shift for large clusters M approximately 
like T c — T m (M) cx 1/M, and at temperatures above T c , 
n m (T) goes to zero inversely proportional to M 2 (for de- 
tails, see Ref. 18). Moreover, the transition is signalled 
by a pronounced maximum in the temperature deriva- 
tive of the energy close to T m , similarly to that of the 
fluctuations in the island size and the energy. All these 
quantities seem to show singular behaviour in the ther- 
modynamic limit. 

In general, a phase diagram in the temperature- 
coverage (kBT/Eb,9) plane may be determined, with a 
single large cluster or, at 9 > 1/4, a stripe of adatoms 
characterising the low-temperature phase. Because of 
the transcription to the Ising model, (3), T c is known 
exactly at 9= 1/2, namely k B T/E b = l/(21n(\/2+ 1)). 
Certainly, the transition temperature tends to decrease 
with decreasing coverage. For instance, we estimated, 
at 9= 0.04, 0.01, and 0.0016 the critical temperatures 
k B T/E b w 0.49, 0.38, and 0.28, respectively. 

Approaching the phase transition from below, the 
largest cluster becomes more and more ramified and may 
dissociate quite easily, while other groups of adatoms may 
coalesce forming a new largest cluster, as one may read- 
ily observe in the simulations. Accordingly, the island 
size n m (T) as well as the energy may fluctuate strongly, 
and good statistics is needed to get reliable equilibrium 
values in the critical region, and to quantify the asymp- 
totics. In particular, we analysed the order parameter 
n m (T), as T — ► T c . Fitting the simulational data to a 
power-law 

n m (T)^(T c -Tf (5) 

we obtained, from extensive simulations at coverages 
0.0016 < 9 < 0.04, S3 = 0.45 ± 0.01, being (if at 
all) only fairly weakly dependent on the coverage. The 
rather large error bar reflects, especially, uncertainties 
in extrapolating the estimates to the thermodynamic 
limit. Actually, we computed an effective exponent 
ef/ = dln(n m )/dln(t), where t =\ T c - T | or | T m - T |, 
leading to upper and lower bounds for the critical expo- 
nent f3, which is approached in the limit L, M — ► co and 
t — > 0. 

The standard description of a cluster in equilibrium 
with a gas is based on the Gibbs-Thomson formula [19]. 
The equilibrium density p of the gas of adatoms in coex- 
istence with a circular island of radius R is then [11,12] 

p{R) = p s exp[j/(R Pi k B T)] (6) 

where 7 is the free energy per unit length of the island 
edge, P i is the density of the island, and p s (T) denotes the 
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density of a gas of adatoms in the presence of a straight 
step. Modified formulae have been discussed [11,14] to 
include, e.g., interactions between the adatoms on the 
terrace. However, no attempt is known to us to extend 
the description into the critical region, where the island 
tends to disappear, p~pi- 

In our simulations, we confirmed that the density 
of adatoms p is, indeed, constant, when exceeding a 
critical distance from the center of the island. More- 
over, at low temperatures and coverages, the Gibbs- 
Thomson formula (6) predicts an approximately loga- 
rithmic dependence of the characteristic temperature T x , 
with n m (T x ) = x, on the coverage 9 

fc B T x oc -l/ln[0(l-a;)] (7) 

fixing the number of adatoms M 2 and varying the surface 
size L 2 . To derive (7), one may use the low-temperature 
relation [12] p s cx exp(— Ead/ksT), where E a d is the en- 
ergy to detach an atom from the step. Indeed, for exam- 
ple, our simulational data at M = 10, 100 < L < 1000, 
and x = 1/2 agree well with (7) [18]. 

IV. DIFFUSION 

The exchange of adatoms between the island and the 
surrounding terrace gas leads to a diffusive motion of the 
equilibrated island. The corresponding diffusion constant 
D follows from the fluctuations in the position of the 
center of mass of the cluster, r CTO (t), 

D=<{v cm {t)-v cm (Q)f> /(At) (8) 

where the brackets denote the equilibrium average. 

From previous studies [1-8], one expects that D de- 
pends on the number of adatoms in the island, N — 
(N max }, at least asymptotically for large values of N, as 

DocN- a (9) 

with the value of a characterising the dominat mecha- 
nism of exchange between the cluster and the gas, dis- 
criminating, e.g., island edge or periphery diffusion, ter- 
race diffusion and evaporation-condensation kinetics (a 
similar classification holds for step fluctuations [20,21]), 
see below. 

To check the diffusive character of the island motion, 
(8), and to determine the characteristic exponent a, we 
performed kinetic Monte Carlo simulations, using the 
bond-counting ansatz, (1), with E n =0, thereby speed- 
ing up the dynamic processes. The time may be mea- 
sured in units of seconds. Invoking the microscopic at- 
tempt frequency v 1 one trial jump per adatom, 1 MCA, 
corresponds then to l/(4i/) seconds [13,16]. We chose 
v = 10 11 Hz. The lattice constant of the square surface 
was set equal to one. To compute the fluctuations in the 
position of the center of mass of the island, (8), we first 
equilibrated the system, before avering over an ensemble 



of initial times as well as an ensemble of (up to 1000) 
realizations. The simulations were done at low tempera- 
tures, well below T c , namely ksT '/Et, = 0.2 and 0.28, for 
islands of sizes N ranging from about 20 to about 700. 
The coverage was fixed, 9 = 0.01. We applied here the 
BKL algorithm in our extensive simulations. 

The positional fluctuations, ((r cm (t) — r cm (0)) 2 ), are 
found to increase, indeed, linearly in time, even at early 
times. The diffusion constant D is then readily obtained 
from linear regression. The resulting size dependent 
D(N) at k B T/E b =0.28 is shown in Fig. 4. Discard- 
ing the smallest island size, the characteristic exponent 
a is estimated to be, on average, a = 1.02 ± 0.03. At 
ksT/Eb^ 0.2, the value of a seems to be consistent with 
a = 1 as well (a = 1.04±0.04), for island sizes N ranging 
from 20 to 700. 

From elementary geometry and energy considerations 
[14] and from a Langevin theory [2], a=l may be argued 
to correspond to an island motion driven by terrace dif- 
fusion, where the adatom emitted by the island diffuses 
as a random walker on the terrace before attaching again 
(or a vacancy may diffuse through the island) . 

In principle, other mechanisms may compete, in partic- 
ular, random detachments and attachments of adatoms 
at the island edge (evaporation-condensation kinetics) 
or diffusion of the atoms along the island edge (periph- 
ery diffusion), leading to a= 1/2 and 3/2, respectively 
[2,14]. In addition, periphery diffusion may be hindered 
by corners and kinks, modifying, possibly, the value of 
a [4,6]. Among othres, details of the shape of the is- 
lands as well as the activation energies are expected to 
determine which mechanism dominates. In general, it is 
reasonable to consider an effective characteristic expo- 
nent a e ff— -dln£>/dlniV, which may depend on island 
size N and temperature T, as observed in experiments 
and simulations. 

In our case, Fig. 4 indicates an increase of a e ff from 
higher values at small island sizes towards a rts 1 at 
TV w 30 — 40. Of course, another crossover at cluster 
sizes exceeding those we studied, N rts 700, cannot be 
excluded. Actually, simulations on a related Ising model, 
applying Kawasaki dynamics, were interpreted as provid- 
ing evidence for a crossover from periphery diffusion to 
evaporation-condensation kinetics, studying somewhat 
smaller cluster sizes, N < 500, and higher coverages, 
with an average exponent a not far from one [14]. In 
contrast to the previous study, we determined D(N) at 
fixed coverage. Obviously, our data do not show clearly 
such a crossover, but it cannot be ruled out. 

It seems interesting to note that a w 1 has been found 
in experiments on Ag(lll), using scanning tunneling mi- 
croscopy [1], but other values of a have been reported for 
different surfaces, reflecting the above mentioned com- 
peting mechanisms. 
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V. SUMMARY 

Thermal stability and diffusion of an adatom island of 
monoatomic height on a square surface have been stud- 
ied, using a standard bond-counting ansatz for the en- 
ergy barriers and performing kinetic Monte Carlo simu- 
lations. 

A phase transition has been identified, with the frac- 
tion of adatoms in the largest cluster n m being the order 
parameter. n m (T) vanishes on approach to the transition 
temperature T c as n m oc| T c -T f with (3 = 0.45±0.10 at 
various small coverages 9. In addition, the energy as well 
as fluctuations in the island size and in the energy exhibit 
singular behaviour at T c . Experimentally, such phenom- 
ena may be observed on a terrace bordered by descending 
steps with large Schwoebel-Ehrlich barriers to allow for 
equilibration of the island with the surrounding gas of 
adatoms on the terrace. 

The diffusion constant D, describing the island mo- 
tion, decreases with the number of adatoms in the is- 
land N like D <x N~ a , with a being close to one for 
an extended range of island sizes at temperatures well 
below T c at 9 = 0.01. a = 1 corresponds to the case 
where the dominant mechanism driving the motion of 
the island is terrace diffusion. An alternative interpre- 
tation invokes a crossover from periphery diffusion to 
evaporation-condensation kinetics . 

Of course, it would be interesting to investigate the 
robustness of our findings against varying, especially, the 
activation energies in a systematic way. 
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Figure Captions 



Fig. 1: Island of monoatomic height in equilibrium with 
the surrounding gas of adatoms. 

Fig. 2: Time dependence of the island size n max (t) for 
M — 50 and L = 250, at various temperatures 
ksT/Eb. The initial configuration of the simula- 
tions is a square cluster with M 2 adatoms. The 
time is measured in trial jumps per adatom (MCA). 

Fig. 3: Temperature dependence of the thermally aver- 
aged island size n m (T) for surfaces with 125 2 sites 
and 25 2 adatoms (open circles: periodic bound- 
ary conditions; full diamonds: reflecting bound- 
aries) as well as 250 2 sites and 50 2 adatoms (open 
squares: periodic boundary conditions). Error bars 
are smaller than symbol sizes. 

Fig. 4: Logarithm of the diffusion constant D vs. loga- 
rithm of number of adatoms in the largest cluster 
N at coverage 9 = 0.01 and temperature kBT/E b 
=0.28. The dashed line corresponds to a = 1. The 
size of the surfaces (L 2 ) ranged from 50 2 to 280 2 in 
the simulations. 



5 




t (in units of 10,000 MCA) 



